install.packages('https://cran.r-project.org/src/contrib/Archive/interflex/interflex_1.0.8.tar.gz', 
                 repos=NULL, 
                 type='source')


library(magrittr)
library(stringr)
library(plyr)
library(tidyverse)
library(interflex)
library(readxl)
library(gridExtra)


rm(list=ls())
home = 'C:/Users/jdt34/Dropbox/VNA_Responsiveness/Short Article/JOP-dataverse/'


permutation.results = paste0(home, 'Data/RI-analyses.Rds') %>%
  readRDS
experimental.results = paste0(home, 'Data/experimental-analyses.Rds') %>%
  readRDS


dv_survey = paste0(home, 'Data/survey-outcomes.xlsx') %>%
  read_xlsx %>%
  mutate(Treatment=factor(x=Treatment, 
                          levels=c('Control',
                                   'Citizen',
                                   'Firm')),
         Missing=as.integer(is.na(Q1)),
         CitizenXProp.Citizen=Citizen*Prop.Citizen, 
         FirmXProp.Firm=Firm*Prop.Firm) %>%
  subset(!is.na(Treatment))
class(dv_survey) = 'data.frame'


tmp = subset(dv_survey, !is.na(Q1)) %>%
  ggplot(aes(x=Prop.Citizen, fill=factor(Citizen))) +
  geom_histogram(color='black', bins=50) +
  scale_fill_manual(values=c('white','black'))
citizen.survey.histogram = layer_data(tmp)[,c('fill','xmin','xmax','ymin','ymax')] %>% 
  mutate(xmax=0.974*(xmax+abs(min(xmin))), xmin=0.974*(xmin+abs(min(xmin))), 
         ymin=ymin/160-1, ymax=ymax/160-1)


tmp = subset(dv_survey, !is.na(Q1)) %>%
  ggplot(aes(x=Prop.Firm, fill=factor(Firm))) +
  geom_histogram(color='black', bins=50) +
  scale_fill_manual(values=c('white','black'))
firm.survey.histogram = layer_data(tmp)[,c('fill','xmin','xmax','ymin','ymax')] %>%
  mutate(xmax=0.974*(xmax+abs(min(xmin))), xmin=0.974*(xmin+abs(min(xmin))),
         ymin=ymin/160-1, ymax=ymax/160-1)
rm(tmp)


rimer = function(result, X, Z, dv) {
  if(Z=='Prop.Citizen') {
    if.Z = 'FirmXProp.Firm'
    Zrange = seq(from=0, 
                 to=0.71, 
                 by=0.01)
    xAxisLabel = 'Moderator: % citizen'
    yAxisLabel = paste('Marginal effect of citizen treatment\non', dv)
    TreatedHistogram = subset(citizen.survey.histogram, 
                              fill=='black')
    ControlHistogram = subset(citizen.survey.histogram, 
                              fill=='white')
  } else if(Z=='Prop.Firm') {
    if.Z = 'CitizenXProp.Citizen'
    Zrange = seq(from=0, 
                 to=0.88, 
                 by=0.01)
    xAxisLabel = 'Moderator: % firm'
    yAxisLabel = paste('Marginal effect of firm treatment\non', dv)
    TreatedHistogram = subset(firm.survey.histogram, 
                              fill=='black')
    ControlHistogram = subset(firm.survey.histogram, 
                              fill=='white')
  }
  me.sigs = NULL
  ri.ests = (subset(permutation.results, 
                    Result==result & term==X)$estimate + 
               (subset(permutation.results, 
                       Result==result & term==paste(X, Z, sep=':'))$estimate %*% t(Zrange)))
  our.est = (subset(experimental.results, 
                    Result==result & term==X)$estimate + 
               subset(experimental.results, 
                      Result==result & term==paste(X, Z, sep=':'))$estimate*Zrange)
  for(z in 1:length(Zrange)) {
    qest = ecdf(ri.ests[,z])
    me.sigs[z] = qest(our.est[z]) %>%
      multiply_by(1e3) %>%
      round %>%
      as.integer %>%
      mapvalues(from=seq(0L, 1000L, 1L),
                to=c(rep('p \u2264 0.01', 6), 
                     rep('p \u2264 0.05', 20), 
                     rep('p > 0.05', 950), 
                     rep('p \u2264 0.05', 20), 
                     rep('p \u2264 0.01', 5)),
                warn_missing=F)
  }
  me.sigs = factor(x=me.sigs, 
                   levels=c('p > 0.05',
                            'p \u2264 0.05',
                            'p \u2264 0.01'))
  Labels = list(expression(p > 0.05), 
                expression(p <= 0.05), 
                expression(p <= 0.01))
  if.D = X
  if.X = Z
  if.Z = c(if.Z, 
           setdiff(subset(experimental.results, 
                          Result==result)$term, 
                   c(if.D, if.X, '(Intercept)','Citizen:Prop.Citizen','Firm:Prop.Firm')))
  kern.est = inter.kernel(Y='Q1', 
                          D=if.D, 
                          X=if.X, 
                          Z=if.Z, 
                          data=subset(dv_survey,
                                      CentNom<0.5), 
                          na.rm=T, 
                          CI=F, 
                          cores=6)$est[,c('X','ME')] %>%
    set_colnames(c('Z','dYdX'))
  ri.ests = as.data.frame(ri.ests) %>%
    mutate(ID=1:1e4) %>%
    gather(key='Z', value='dYdX', -ID) %>% 
    mutate(Z=mapvalues(x=Z, 
                       from=paste0('V', 1:length(Zrange)), 
                       to=Zrange),
           Z=as.numeric(Z))
  data.frame(Significance=me.sigs, 
             Z=Zrange, 
             dYdX=our.est) %>%
    ggplot(aes(x=Z, y=dYdX, fill=Significance)) +
    geom_line(data=ri.ests, 
              aes(x=Z, y=dYdX, group=ID), 
              color='black', alpha=0.01, inherit.aes=F) +
    geom_line(data=kern.est, 
              aes(x=Z, y=dYdX), 
              color='#e63e40', linetype=1, alpha=1, size=0.5, inherit.aes=F) + 
    geom_hline(yintercept=0, linetype=2, color='white') +
    geom_point(size=2, shape=21) +
    geom_rect(data=TreatedHistogram, 
              aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), 
              fill='black', color='black', size=0.25, inherit.aes=F) +
    geom_rect(data=ControlHistogram, 
              aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), 
              fill='white', color='black', size=0.25, inherit.aes=F) +
    scale_x_continuous(breaks=seq(from=0, 
                                  to=1, 
                                  by=0.2), 
                       labels=paste0(str_pad(string=seq(from=0, 
                                                        to=100, 
                                                        by=20), 
                                             width=3, 
                                             side='left', 
                                             pad=' '), 
                                     '%')) +
    scale_fill_brewer(labels=Labels, type='seq', palette='YlGnBu', direction=1, drop=F) +
    guides(fill=guide_legend(override.aes=list(size=2))) +
    labs(x=xAxisLabel, y=yAxisLabel) +
    coord_fixed(ratio=0.25, xlim=c(0,1), ylim=c(-1,1), expand=F) +
    theme_minimal() +
    theme(axis.text=element_text(size=8, color='black'),
          axis.title=element_text(size=9, color='black'),
          strip.text=element_text(size=9),
          legend.title=element_text(size=9),
          legend.text=element_text(size=8),
          legend.position='bottom',
          plot.margin=unit(c(0.2, 0.2, 0, 0.1), 'in'),
          panel.border=element_rect(color='black', fill=NA))   
}


fig_9_1 = rimer(result='tab2col7', 
                X='Citizen', 
                Z='Prop.Citizen', 
                dv='debate preparation')


ggsave(filename='figure-A09-1.png', 
       plot=fig_9_1, 
       path=paste0(home, 'Figures/'), 
       width=4.5, 
       height=3, 
       units='in')
ggsave(filename='figure-A09-1.eps', 
       plot=fig_9_1, 
       path=paste0(home, 'Figures/'), 
       width=4.5, 
       height=3, 
       units='in', 
       device=cairo_ps)

